*******************************************************************************************
******** Summary Stats
*******************************************************************************************
clear
eststo clear
set more off
use hourlypanel.dta

keep if year < 2013
g runvarabs = abs(run_var_end)
keep if runvarabs < 35
g midnight = hour == 0

bysort date_m: egen dc_total_inc = sum(psa_num_inc)
bysort date_m eleven: egen dc_total_inc_11 = sum(psa_num_inc) 
bysort date_m midnight: egen dc_total_inc_12b = sum(psa_num_inc) 

bysort psa date_m: egen psa_day_inc = sum(psa_num_inc) 
bysort psa date_m eleven: egen psa_day_inc_11 = sum(psa_num_inc)
bysort psa date_m midnight: egen psa_day_inc_12b = sum(psa_num_inc)


foreach c in inc {
replace psa_day_`c'_12b = -10 if midnight == 0
bysort psa date_m: egen psa_day_`c'_12 = max(psa_day_`c'_12b)

replace dc_total_`c'_12b = -10 if midnight ==0 
bysort date_m: egen dc_total_`c'_12 = max(dc_total_`c'_12b)
}

keep if eleven == 1 

drop psa_day_inc_12b dc_total_inc_12b 

estpost sum psa_day*
#delimit ;
estout using sum_psa_sst.tex, cells("count(fmt(0)) mean(fmt(3)) sd(fmt(3)) min(fmt(0)) max(fmt(0))") nonumber style(tex) replace mlabels(none) collabels(none)
	varlabels(psa_day_inc_12 "Daily PSA SST-detected incidents, midnight-1am" psa_day_inc "Daily PSA SST-detected incidents" psa_day_inc_11 "Daily PSA SST-detected incidents, 11pm-midnight" )
	order(psa_day_inc psa_day_inc_11 psa_day_inc_12)
	postfoot("\hline" "\multicolumn{6}{l}{\textbf{Crime at the PSA-level (geographic unit of analysis)}}\\")
	
	;
#delimit cr	


keep if psa == 701
estpost sum dc*
#delimit ;
estout using sum_dc_sst.tex, cells("count(fmt(0)) mean(fmt(3)) sd(fmt(3)) min(fmt(0)) max(fmt(0))") nonumber style(tex) replace mlabels(none) collabels(none)
	varlabels(dc_total_inc_12 "Daily DC SST-detected incidents, midnight-1am" dc_total_inc "Daily DC SST-detected incidents" dc_total_inc_11 "Daily DC SST-detected incidents, 11pm-midnight" )
	prehead("\begin{table}[htbp]\centering" "\footnotesize" "\vspace{.3cm}" "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" "\caption{Summary Statistics} \label{tab:sumstat}" "\begin{tabular}{l*{5}{r}}" "\hline""\hline") 
	posthead("& N & Mean	& SD & Min	&    Max   \\""\hline"" \multicolumn{6}{l}{\textbf{Gunshots in Washington, DC}}\\")
	postfoot("\hline" "\multicolumn{6}{l}{\textbf{Gunshots at the PSA-level (geographic unit of analysis)}}\\")
	order(dc_total_inc dc_total_inc_11 dc_total_inc_12 )
	;
#delimit cr	

**
clear
eststo clear
set more off
use hourlypanel.dta

keep if year > 2010 & year < 2013
g runvarabs = abs(run_var_end)
keep if runvarabs < 35

rename num_calls psa_num_calls
rename num_calls_mpd psa_num_calls_mpd
rename num_calls_gunshot psa_num_calls_gun


foreach c in rep hom gun vio calls calls_mpd calls_gun{
bysort psa date_m: egen psa_day_`c' = sum(psa_num_`c')
bysort psa date_m eleven: egen psa_day_`c'_11 = sum(psa_num_`c')
}

keep if eleven == 1
estpost sum psa_day*
#delimit ;
estout using sum_psa_crime.tex, cells("count(fmt(0)) mean(fmt(3)) sd(fmt(3)) min(fmt(0)) max(fmt(0))") nonumber style(tex) replace mlabels(none) collabels(none)
	varlabels(psa_day_rep "Daily PSA MPD reported crimes" psa_day_rep_11 "Daily PSA MPD reported crimes, 11pm-midnight" psa_day_hom "Daily PSA MPD reported homicides" psa_day_hom_11 "Daily PSA MPD reported homicides, 11pm-midnight" psa_day_gun "Daily PSA MPD reported gun crimes" psa_day_gun_11 "Daily PSA MPD reported gun crimes, 11pm-midnight" 
	psa_day_vio "Daily PSA MPD reported violent crimes" psa_day_vio_11 "Daily PSA MPD reported violent crimes, 11pm-midnight" psa_day_calls "Daily PSA 911 calls" psa_day_calls_11 "Daily PSA 911 calls, 11pm-midnight" 
	psa_day_calls_mpd "Daily PSA 911 calls for police" psa_day_calls_mpd_11 "Daily PSA 911 calls for police, 11pm-midnight" psa_day_calls_gun "Daily PSA 911 calls reporting gunshot" psa_day_calls_gun_11 "Daily PSA 911 calls reporting gunshot, 11pm-midnight")
postfoot("\hline \hline" "\end{tabular}" )
	;
#delimit cr		


******************************************************************************************
******** Figures
******************************************************************************************

clear
eststo clear
set more off
use hourlypanel.dta


keep if year < 2013 & year > 2008

*** Section 1 - Appendix figures centered on July 4 
bys year doy: egen total_shots = sum(psa_num_inc)
bys doy: egen mean_dc_shots = mean(total_shots)


* Figure A1
#delimit ;
scatter mean_dc_shots doy, msymbol(Oh)
	graphregion(color(white) lc(white))
	yti("Mean daily gunshots")
	xti("Day of year", margin(t=2 b=2))
	xlab(1 "Jan 1" 60 "Mar 1" 182 "July 1" 244 "Sept 1" 121 "May 1"  305 "Nov 1" 365 "Dec 31")
	legend(off)
;
graph export july4_allyear.png, replace;
#delimit cr

*Figure A2
#delimit ;
scatter mean_dc_shots doy if (abs(run_var_start) < 60) & mean_dc_shots < 65, xline(181.5) msymbol(Oh)
	graphregion(color(white) lc(white))
	yti("Mean daily gunshots")
	xti("Day of year", margin(t=2 b=2))
	xlab(182 "July 1" 152 "June 1" 213 "Aug 1" 244 "Sept 1" 121 "May 1")
	legend(off)
;
graph export july4_zoomed.png, replace;
#delimit cr


***********************************************************************************************
*** Main Tables
***********************************************************************************************

clear
eststo clear
set more off
use hourlypanel.dta

*** Section 2 - Main figure and Appendix Figures 3, 4 and 5
keep if hour == 23 | hour == 0

*weekday indicator (0 is sunday, 1 is monday)
g wkday = 0
	replace wkday = 1 if dow < 5

*two way interactions	
g curfew_wkday = wkday * curfew_early
g curfew_11 = eleven * curfew_early
g wkday_11 = wkday * eleven

*three way interaction (for Triple Diff)
g dddtreat = wkday*curfew_early*eleven
	
*run var and interaction 
g run_var = run_var_end
g rv_curfew = run_var*curfew_early

g run_var_both = run_var_start
		replace run_var_both = run_var_end if end == 1
g rv_both_curfew = run_var_both*curfew_early
g rv_s = run_var_both*start
g rv_both_curfew_s = run_var_both*curfew_early*start

foreach v in run_var rv_curfew {
g `v'_2 = `v'^2
g `v'_3 = `v'^3
}

*making running var interactions
foreach var in wkday curfew_wkday eleven curfew_11 dddtreat wkday_11{
g rv_`var' = run_var * `var'
g rv_`var'_2 = (run_var * `var')^2
g rv_`var'_3 = (run_var * `var')^3

g rv_`var'_both = run_var_both * `var'
}
foreach var in wkday curfew_wkday eleven curfew_11 dddtreat wkday_11{
g rv_`var'_both_s = run_var_both * `var'* start
}
*dow: zero is sunday
forvalues d = 0/6{
g dow_`d' = 0
	replace dow_`d' = 1 if dow == `d'
}

g july17_dum = 0
	replace july17_dum = 1 if doy > 181 & doy < 189

egen wkbin = cut(run_var), at(-35(7)35)	
*drop if abs(run_var) > 33


global fe_weather prcp tmax tmin school i.psa i.year dow_1 dow_2 dow_3 dow_4
*dow dummies hardcoded because i.dow can lead to dropping wkday

global dddvars dddtreat curfew_early eleven wkday curfew_wkday curfew_11 wkday_11 run_var rv_curfew rv_wkday rv_eleven rv_curfew_wkday rv_curfew_11 rv_wkday_11 rv_dddtreat 
global dddvars_2 dddtreat curfew_early eleven wkday curfew_wkday curfew_11 wkday_11 run_var rv_curfew rv_wkday rv_eleven rv_curfew_wkday rv_curfew_11 rv_wkday_11 rv_dddtreat run_var_2 rv_curfew_2 rv_wkday_2 rv_eleven_2 rv_curfew_wkday_2 rv_curfew_11_2 rv_wkday_11_2 rv_dddtreat_2 
global dddvars_3 dddtreat curfew_early eleven wkday curfew_wkday curfew_11 wkday_11 run_var rv_curfew rv_wkday rv_eleven rv_curfew_wkday rv_curfew_11 rv_wkday_11 rv_dddtreat run_var_2 rv_curfew_2 rv_wkday_2 rv_eleven_2 rv_curfew_wkday_2 rv_curfew_11_2 rv_wkday_11_2 rv_dddtreat_2 run_var_3 rv_curfew_3 rv_wkday_3 rv_eleven_3 rv_curfew_wkday_3 rv_curfew_11_3 rv_wkday_11_3 rv_dddtreat_3  
global dddvars_start dddtreat curfew_early start eleven wkday curfew_wkday curfew_11 wkday_11 run_var_both rv_s rv_both_curfew rv_both_curfew_s rv_wkday_both rv_wkday_both_s rv_eleven_both rv_eleven_both_s rv_curfew_wkday_both rv_curfew_wkday_both_s rv_curfew_11_both rv_curfew_11_both_s rv_wkday_11_both rv_wkday_11_both_s rv_dddtreat_both rv_dddtreat_both_s
global dirdvars_end curfew_wkday wkday curfew_early run_var rv_curfew rv_wkday rv_curfew_wkday


eststo r1: xi: reg psa_num_inc $dddvars $fe_weather if (run_var_end > -34 & run_var_end < 34) & year <2013, vce(cluster run_var_end)
eststo r2: xi: reg psa_num_inc $dddvars_start $fe_weather if (run_var_both > -32 & run_var_both < 32) & july17_dum == 0 & year <2013, vce(cluster doy) /*Cannot use 34 days bandwidth because only 62 days between cutoffs, 31 is max */
eststo r3: xi: reg psa_num_inc $dddvars_2 $fe_weather if (run_var_end > -34 & run_var_end < 34) & year <2013, vce(cluster run_var_end)
eststo r4: xi: reg psa_num_inc $dddvars_3 $fe_weather if (run_var_end > -34 & run_var_end < 34) & year <2013, vce(cluster run_var_end)

#delimit ;
estout r1 r2 r3 r4 using table2.tex, replace
	posthead( "\hline" "\\""\\[-.5cm]")
	prehead("\begin{table}[htbp]\centering" "\footnotesize" "\vspace{.3cm}" "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" "\caption{Effect of the 11pm curfew on gunfire incidents}" "\label{tab:curfew_sst}" "\begin{tabular}{l*{@M}{c}}" "\hline" "\hline" "&Main&Include July&Quadratic&Cubic\\" "&specification&curfew change&time trend&time trend\\" "&(1)&(2)&(3)&(4)\\")
	prefoot("\\""\\[-.5cm]" "\hline" "\\""\\[-.5cm]") collabels(none) mlabels(none)
	postfoot( "\\""\\[-.5cm]" "\hline" "\hline"  "\end{tabular}" ) keep(dddtreat) order(dddtreat)
	c(b(star fmt(3)) se(par fmt(3))) starlevels(* 0.1 ** 0.05 *** 0.01) varlabels(dddtreat "Early curfew * Weekday * 11pm hour")
	stats(N, fmt(0) label("Observations")) style(tex)
	;
#delimit cr

#delimit ;
estout r1 using ddd_allvars.tex, replace
posthead( "\hline" "\\""\\[-.5cm]")
	prehead("\begin{table}[htbp]\centering" "\footnotesize" "\vspace{.3cm}" "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" "\caption{D-i-D-i-RD}" "\begin{tabular}{l*{@M}{c}}" "\hline""\hline") 
	prefoot("\\""\\[-.5cm]" "\hline" "\\""\\[-.5cm]") mlabels("1") collabels(none) 
	postfoot( "\\""\\[-.5cm]" "\hline" "\hline"  "\end{tabular}") 
	keep(dddtreat curfew_early eleven wkday curfew_wkday curfew_11 wkday_11 run_var rv_curfew rv_wkday rv_eleven rv_curfew_wkday rv_curfew_11 rv_wkday_11 rv_dddtreat prcp tmax tmin school) 
	order(dddtreat wkday curfew_early curfew_wkday eleven curfew_11 wkday_11 )
	c(b(star fmt(4)) se(par fmt(4))) starlevels(* 0.1 ** 0.05 *** 0.01) 
	varlabels(dddtreat "Early curfew * Weekday * 11pm hour" eleven "11pm hour" curfew_11 "Early curfew * 11pm hour" wkday_11 "Weekday * 11pm hour" curfew_wkday "Early curfew * Weekday" wkday "Weekday" curfew_early "Early curfew"
	run_var "Time" rv_curfew "Time * Early curfew" rv_wkday "Time * Weekday" rv_eleven "Time * 11pm hour" rv_curfew_wkday "Time * Early curfew * Weekday" rv_curfew_11 "Time * Early curfew * 11pm hour" 
	rv_wkday_11 "Time * Weekday * 11pm hour" rv_dddtreat "Time * Early curfew * Weekday * 11pm hour" prcp "Precipitation (mm)" tmax "Max temp (Celsius)" tmin "Min temp (Celsius)" school "School in session")
	stats(N, fmt(0) label("Observations")) style(tex)
;
#delimit cr

eststo clear

foreach c in psa_num_inc psa_num_rep psa_num_hom psa_num_gun psa_num_vio num_calls num_calls_mpd num_calls_gunshot{
eststo: xi: reg `c' $dirdvars_end $fe_weather if (run_var_end > -34 & run_var_end < 34) & eleven == 1  & year > 2010 & year <2013, vce(cluster run_var_end) 
}


#delimit ;
estout using table3.tex, replace 
	posthead( "\hline" "\\""\\[-.5cm]")
	prehead("\begin{table}[htbp]\centering" "\footnotesize" "\vspace{.3cm}" "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" "\caption{Effect of the 11pm curfew on other measures of gun violence} \label{tab:compare_sst}" "\begin{tabular}{l*{8}{c}}"
		"\hline" "\hline"
		"&\multicolumn{1}{c}{}&\multicolumn{4}{c}{\underline{MPD Reported Crime}}&\multicolumn{3}{c}{\underline{911 Calls for Service}}\\"
		"&\multicolumn{1}{c}{ShotSpotter}&\multicolumn{1}{c}{}&\multicolumn{1}{c}{}&\multicolumn{1}{c}{Gun-}&\multicolumn{1}{c}{}&\multicolumn{1}{c}{}&\multicolumn{1}{c}{For}&\multicolumn{1}{c}{To Report}\\"
		"&\multicolumn{1}{c}{Incidents}&\multicolumn{1}{c}{All}&\multicolumn{1}{c}{Homicides}&\multicolumn{1}{c}{Involved}&\multicolumn{1}{c}{Violent}&\multicolumn{1}{c}{All}&\multicolumn{1}{c}{Police}&\multicolumn{1}{c}{Gunshots}\\"
		"           &         (1)   &         (2)   &         (3)   &         (4)   &         (5)   &         (6)   &         (7)   &         (8)   \\")
	
	prefoot("\\""\\[-.5cm]" "\hline" "\\""\\[-.5cm]") mlabels(none) collabels(none) 
	postfoot("Late curfew mean&0.017&0.100&0.000&0.012&0.04&0.885&0.721&0.013\\" "\\""\\[-.5cm]" "\hline" "\hline"  "\end{tabular}") 
	c(b(star fmt(4)) se(par fmt(4))) starlevels(* 0.1 ** 0.05 *** 0.01)
	style(tex) stats(N, fmt(0) label("Observations"))
	keep(curfew_wkday) 
	varlabels(curfew_wkday "Early curfew * Weekday")
;
#delimit cr

eststo clear
eststo: xi: reg psa_num_inc rv_eleven run_var eleven $fe_weather if dow < 5 & year < 2013 & run_var_end < 0 & run_var_end > -34, vce(cluster run_var_end) 
eststo: xi: reg psa_num_inc rv_wkday run_var wkday $fe_weather if eleven == 1 & year < 2013 & run_var_end < 0 & run_var_end > -34, vce(cluster run_var_end) 

#delimit ;
estout using trendstest_sst.tex, replace
	posthead( "\hline" "\\""\\[-.5cm]")
	prehead("\begin{table}[htbp]\centering" "\footnotesize" "\vspace{.3cm}" "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" "\caption{Tests for diverging trends in ShotSpotter data}" "\begin{tabular}{l*{@M}{c}}" "\hline""\hline") 
	prefoot("\\""\\[-.5cm]" "\hline" "\\""\\[-.5cm]") mlabels("11pm vs. Midnight" "Weekday vs. Weekend" ) collabels(none) 
	postfoot( "\\""\\[-.5cm]" "\hline" "\hline"  "\end{tabular}" ) keep(rv_wkday rv_eleven) order(rv_eleven rv_wkday)
	c(b(star fmt(4)) se(par fmt(4))) starlevels(* 0.1 ** 0.05 *** 0.01) varlabels(rv_eleven "11pm * Time" rv_wkday "Weekday * Time" )
	stats(N, fmt(0) label("Observations")) style(tex)
;
#delimit cr

eststo clear
foreach c in psa_num_inc psa_num_rep num_calls {
eststo: xi: reg `c' rv_wkday run_var wkday $fe_weather if eleven ==1 & run_var_end < 0 & run_var_end > -34 & year > 2010,  vce(cluster run_var_end)
}

#delimit ;
estout using trendstest_compwkday.tex, replace
	posthead( "\hline" "\\""\\[-.5cm]" "\textbf{Weekday vs. Weekend}&&&\\")
	prehead("\begin{table}[htbp]\centering" "\footnotesize" "\vspace{.3cm}" "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}" "\caption{Tests for diverging trends in reported crime and 911 call data}" "\begin{tabular}{l*{@M}{c}}" "\hline""\hline") 
	mlabels("SST" "Reported crime" "911 calls") collabels(none) 
	postfoot( "\hline" ) keep(rv_wkday) order(rv_wkday)
	c(b(star fmt(4)) se(par fmt(4))) starlevels(* 0.1 ** 0.05 *** 0.01) varlabels(rv_wkday "Weekday * Time" )
	stats(N, fmt(0) label("Observations")) style(tex)
;
#delimit cr
eststo clear


foreach c in psa_num_inc psa_num_rep num_calls {
eststo: xi: reg `c' rv_eleven run_var eleven $fe_weather if dow < 5 & run_var_end < 0 & run_var_end > -34 & year > 2010 , vce(cluster run_var_end)
}


#delimit ;
estout using trendstest_comp11.tex, replace
	posthead("\textbf{11pm vs. Midnight}&&&\\")
	mlabels(none)  collabels(none) 
	postfoot( "\\""\\[-.5cm]" "\hline" "\hline"  "\end{tabular}") keep(rv_eleven) order(rv_eleven)
	c(b(star fmt(4)) se(par fmt(4))) starlevels(* 0.1 ** 0.05 *** 0.01) varlabels(rv_eleven "11pm*running var")
	stats(N, fmt(0) label("Observations")) style(tex)
;
#delimit cr


keep if year <2013
drop if abs(run_var_end) > 33

global fe_weather prcp tmax tmin school i.psa i.year i.dow
xi: reg psa_num_inc $fe_weather run_var, vce(cluster run_var)  
predict resid_2, r

g r_wkday_11 = resid_2 * wkday *eleven
g r_wkend_11 = resid_2 * (1-wkday) *eleven
g r_wkday_12 = resid_2 * wkday * (1-eleven)
g r_wkend_12 = resid_2 * (1-wkday) *(1-eleven)

drop wkday_11
g wkday_11 = wkday *eleven * psa_num_inc
	replace wkday_11 = . if wkday == 0 | eleven == 0
g wkend_11 = (1-wkday) *eleven * psa_num_inc
	replace wkend_11 = . if wkday == 1 | eleven == 0
g wkday_12 = wkday * (1-eleven) * psa_num_inc
	replace wkday_12 = . if wkday == 0 | eleven == 1
g wkend_12 = (1-wkday) *(1-eleven) * psa_num_inc
	replace wkend_12 = . if wkday == 1 | eleven == 1

foreach r in r_wkday_11 r_wkend_11 r_wkday_12 r_wkend_12{
bys run_var: egen test1 = sum(`r') 
bys wkbin: egen test1b = sum(`r') 
replace `r' = . if `r' ==0
bys run_var: egen test2 = count(`r') 
bys wkbin: egen test2b = count(`r')
g `r'_mean = test1/test2
g `r'_mean_wk = test1b/test2b
drop test1 test2 test1b test2b
}
foreach r in wkday_11 wkend_11 wkday_12 wkend_12{
bys run_var: egen test1 = sum(`r') 
bys run_var: egen test2 = count(`r') 
g `r'_mean = test1/test2
drop test1 test2
}

*keep run_var r_wkday_11_mean r_wkend_11_mean r_wkday_12_mean r_wkend_12_mean r_wkday_11_mean_wk r_wkend_11_mean_wk r_wkday_12_mean_wk r_wkend_12_mean_wk wkbin wkday_11_mean wkend_11_mean wkday_12_mean wkend_12_mean
duplicates drop

g ddd = (r_wkday_11_mean - r_wkend_11_mean) - (r_wkday_12_mean - r_wkend_12_mean)

g new_res_wkday = r_wkday_11_mean - r_wkday_12_mean

g new_res_11 = r_wkday_11_mean - r_wkend_11_mean 


*Figure 1
#delimit ;
twoway (scatter ddd run_var_end if abs(run_var_end) < 34, xline(-.5) msymbol(Oh))
(lfit ddd run_var_end if run_var_end < 0)
(lfit ddd run_var_end if run_var_end > -1),
	xlab( 0 "Sept 1" 15 "Sept 15" 31 "Oct 1" -17 "Aug 15" -31 "Aug 1")
	t1("late curfew                                         early curfew", margin(t=4))
	xti("Day of year", margin(t=2 b=2))
	yti("Mean residuals")
	legend(off)
	graphregion(color(white) lc(white))
;
graph export ddd_end.png, replace;
#delimit cr


*Figure A3
#delimit ;
twoway (scatter new_res_11 run_var_end if abs(run_var_end) < 34, xline(-.5) msymbol(Oh))
(lfit new_res_11 run_var_end if run_var_end < 0)
(lfit new_res_11 run_var_end if run_var_end > -1),
	xlab( 0 "Sept 1" 15 "Sept 15" 31 "Oct 1" -17 "Aug 15" -31 "Aug 1")
	t1("late curfew                                         early curfew", margin(t=4))
	xti("Day of year", margin(t=2 b=2))
	yti("Mean residuals")
	legend(off)
	graphregion(color(white) lc(white))
;
graph export new_res_11.png, replace;
#delimit cr

*Figure A4
#delimit ;
twoway (scatter new_res_wkday run_var_end if abs(run_var_end) < 34, xline(-.5) msymbol(Oh))
(lfit new_res_wkday run_var_end if run_var_end < 0)
(lfit new_res_wkday run_var_end if run_var_end > -1),
	xlab( 0 "Sept 1" 15 "Sept 15" 31 "Oct 1" -17 "Aug 15" -31 "Aug 1")
	t1("late curfew                                         early curfew", margin(t=4))
	xti("Day of year", margin(t=2 b=2))
	yti("Mean residuals")
	legend(off)
	graphregion(color(white) lc(white))
;
graph export new_res_wkday.png, replace;
#delimit cr


*Figure A5 (top)
*
twoway (scatter r_wkday_11_mean run_var if run_var <0) (scatter r_wkend_11_mean run_var if run_var <0, msymbol(oh)) (lpoly resid_2 run_var if run_var <0 & eleven == 1 & wkday == 1 ) /*
	*/ (lpoly resid_2 run_var if run_var <0 & eleven == 1 & wkday == 0, lp(dash)), /* 
	*/ legend(lab(1 "11pm weekday") lab(2 "11pm weekend") lab(3 "11pm weekday") lab(4 "11pm weekend") order(1 2 3 4)) graphregion(color(white) lc(white)) /*
	*/ xlab( 0 "Sept 1" -17 "Aug 15" -31 "Aug 1") ysc(r(-.06 .16)) ylab(-.05 0 .05 .1 .15) yti("Mean residuals") xti("Day of year", margin(t=2 b=2))
graph export resid_dird_lpoly_pre.png, replace
*
*Figure A5 (bottom)
*
twoway (scatter r_wkday_11_mean run_var if run_var <0) (scatter r_wkday_12_mean run_var if run_var <0, msymbol(oh)) (lpoly resid_2 run_var if run_var <0 & eleven == 1 & wkday == 1 ) /*
	*/ (lpoly resid_2 run_var if run_var <0 & eleven == 0 & wkday == 1, lp(dash)), ysc(r(-.06 .16))/* 
	*/ legend(lab(1 "11pm weekday") lab(2 "12am weekday") lab(3 "11pm weekday") lab(4 "12am weekday") order(1 2 3 4))  /*
	*/ graphregion(color(white) lc(white)) ylab(-.05 0 .05 .1 .15) xlab( 0 "Sept 1" -17 "Aug 15" -31 "Aug 1") yti("Mean residuals")	xti("Day of year", margin(t=2 b=2))
graph export resid_ddd_lpoly_pre.png, replace
*

******************************************************************************************
*Section 3 - Lags and Leads figure

clear
eststo clear
set more off
use hourlypanel.dta

keep if year < 2013
keep if hour == 23 | hour == 0

g run_var = run_var_end
g wkday = 0
	replace wkday = 1 if dow < 5
		
*two way interactions	
g curfew_wkday = wkday * curfew_early
g curfew_11 = eleven * curfew_early
g wkday_11 = wkday * eleven

*three way interaction (for Triple Diff)
g dddtreat = wkday*curfew_early*eleven

*run var and interaction 
g rv_curfew = run_var*curfew_early

foreach v in run_var rv_curfew{
g `v'_2 = `v'^2
}

*making running var interactions
foreach var in wkday curfew_wkday eleven curfew_11 dddtreat wkday_11{
g rv_`var' = run_var * `var'
g rv_`var'_2 = (run_var * `var')^2
}

forvalues i = 1/14{
g lead`i' = 0
	replace lead`i' = 1 if wkday == 1 & run_var == - `i' & eleven == 1 
}

forvalues i = 0/14{
g lag`i' = 0
	replace lag`i' = 1 if wkday == 1 & run_var == `i' & eleven == 1 
}

g lag15 = 0
	replace lag15 = 1 if run_var>14 & wkday == 1 & eleven == 1 
****** making the fixed effects: rdrobust can't use the automatic ones


*dow: zero is sunday
forvalues d = 0/6{
g dow_`d' = 0
	replace dow_`d' = 1 if dow == `d'
}

global dddvars_end lag* lead* curfew_early eleven wkday curfew_wkday curfew_11 wkday_11 run_var rv_curfew rv_wkday rv_eleven rv_curfew_wkday rv_curfew_11 rv_wkday_11 rv_dddtreat 
global dirdvars_end lag* lead* wkday curfew_early run_var rv_curfew rv_wkday rv_curfew_wkday

global fedow dow_1 dow_2 dow_3 dow_4 dow_5
global fe_weather prcp tmax tmin school i.psa i.year $fedow


eststo clear
set more off

eststo: xi: reg psa_num_inc $dddvars_end $fe_weather if (abs(run_var) < 34), vce(cluster doy)

forvalues i = 1/14{
g b_lead_`i' = _b[lead`i']
g se_lead_`i' = _se[lead`i']
}
forvalues i = 0/15{
g b_lag_`i' = _b[lag`i']
g se_lag_`i' = _se[lag`i']
}

keep se* b*

duplicates drop

g id = 1

reshape long b_lag_ se_lag_ b_lead_ se_lead_ , i(id) j(band, string)

destring band, replace force

drop id 
save ddd_leads, replace
foreach b in lag lead{
g ul_`b' = b_`b'_ + (1.96*se_`b'_)
g ll_`b' = b_`b'_ - (1.96*se_`b'_)
g ul_`b'_90 = b_`b'_ + (1.65*se_`b'_)
g ll_`b'_90 = b_`b'_ - (1.65*se_`b'_)
}

g band_lead = band*-1

*Figure 2
#delimit ;
twoway rcap ll_lead ul_lead band_lead, lc(red)|| rcap ll_lead_90 ul_lead_90 band_lead, lc(black) || 
	rcap ll_lag ul_lag band, lc(red)|| rcap ll_lag_90 ul_lag_90 band, lc(black) || 
	scatter  b_lag_ band, mc(green) || scatter  b_lead_ band_lead, mc(green) yline(0, lcolor(black)) xline(-.5, lcolor(red))
	legend(order(1 2 5) label(1 "95% CI") label( 2 "90% CI") label( 5 "Coefficient")) 
	xtitle("Days from curfew change", margin(t=2 b=2)) ytitle("Estimated effect", margin(l=2 r=2))
	graphregion(color(white) lc(white))
	;
graph export ddd_lags_leads15.png, replace;
#delimit cr


**************************************************************************************************
**************************************************************************************************
**************************************************************************************************
** Section 4 - makes figure 3 (effect on other hours of the day)

clear
eststo clear
set more off
use hourlypanel.dta

keep if year < 2013

*weekday indicator
g wkday = 0
	replace wkday = 1 if dow < 5
	
*two way interactions	
g curfew_wkday = wkday * curfew_early
g curfew_11 = eleven * curfew_early
g wkday_11 = wkday * eleven

*run var and interaction 
g run_var = run_var_end 
g rv_curfew = run_var*curfew_early

foreach v in run_var rv_curfew{
g `v'_2 = `v'^2
}

*making running var interactions
foreach var in wkday curfew_wkday eleven curfew_11 wkday_11{
g rv_`var' = run_var * `var'
g rv_`var'_2 = (run_var * `var')^2
}

*dow: zero is sunday
forvalues d = 0/6{
g dow_`d' = 0
	replace dow_`d' = 1 if dow == `d'
}

*creating globals for sets of fes - omitting the first one for each group
global fepsa psa_501 psa_601 psa_701 psa_302 psa_502 psa_602 psa_702 psa_303 psa_503 psa_603 psa_703 psa_304 psa_504 psa_604 psa_704 psa_305 psa_505 psa_605 psa_705 psa_306 psa_506 psa_606 psa_706 psa_307 psa_507 psa_607 psa_707 psa_308 psa_608 psa_708
global feyear year_2007 year_2008 year_2009 year_2010 year_2011 year_2012
global fedow dow_1 dow_2 dow_3 dow_4 dow_5 dow_6

global dirdvars_end curfew_wkday wkday curfew_early run_var rv_curfew rv_wkday rv_curfew_wkday

global fe_weather prcp tmax tmin school i.year dow_1 dow_2 dow_3 dow_4

eststo clear

estimates clear

forvalues h = 0/23{
xi: areg psa_num_inc $dirdvars_end $fe_weather if hour == `h' & abs(run_var_end) <34, vce(cluster doy) absorb(psa)
g beta_e_dird_`h' = _b[curfew_wkday]
g se_e_dird_`h' = _se[curfew_wkday]
}

keep beta* se*
duplicates drop

g id = 1

#delimit ;
reshape long beta_e_dird_ se_e_dird_ 
beta_e_uni_ se_e_uni_ beta_e_tri_ se_e_tri_ 
, i(id) j(hour)
;
#delimit cr


g ul_e_dird = beta_e_dird_ + (1.96*se_e_dird_ )
g ll_e_dird = beta_e_dird_ - (1.96*se_e_dird_ )
g ul90_e_dird = beta_e_dird_ + (1.65*se_e_dird_ )
g ll90_e_dird = beta_e_dird_ - (1.65*se_e_dird_ )

g hour_fig = hour
replace hour_fig = hour + 24 if hour < 6

*Figure 3
#delimit ;
twoway rcap ll_e_dird ul_e_dird hour_fig if hour < 23, legend(label(1 "95% CI")) || rcap ll90_e_dird ul90_e_dird hour_fig if hour < 23, legend(label(2 "90% CI")) || 
	rcap ll_e_dird ul_e_dird hour_fig if hour == 23, legend(label(3 "11pm (switching hour)")) lc(red) || rcap ll90_e_dird ul90_e_dird hour_fig if hour == 23, lc(red) ||
	scatter beta_e_dird_ hour_fig, yline(0, lcolor(black))  legend( order(1 2 3 5) label(5 "Coefficient")) xtitle("Hour of day", margin(t=2 b=2)) ytitle("Estimated effect", margin(l=2 r=2))
		graphregion(color(white) lc(white))  xlab(6 "6 am" 12 "noon" 18 "6 pm" 24 "midnight")
	;
graph export hours_e_dird.png, replace;
#delimit cr


